d02kef
d02kef
© Numerical Algorithms Group, 2002.
Purpose
D02KEF 2nd order Sturm-Liouville problem, regular/singular system,
finite/infinite range, eigenvalue and eigenfunction, user-specified
break-points
Synopsis
[elam,hmax,delam,match,maxit,ifail] = d02kef(xpoint,coeffn,bdyval,k,tol,...
elam,monit,report<,hmax,delam,match,maxfun,maxit,ifail>)
Description
~~~~~~~~
D02KEF finds a specified eigenvalue (lambda) of a Sturm-
Liouville system defined by a self-adjoint differential equation
of the second-order
(p(x)y')'+q(x;(lambda))y=0, a<x<b
together with the appropriate boundary conditions at the two
(finite or infinite) end-points a and b. The functions p and q,
which are real-valued, must be defined by a subroutine COEFFN.
The boundary conditions must be defined by a subroutine BDYVAL,
and, in the case of a singularity at a or b, take the form of an
asymptotic formula for the solution near the relevant end-point.
~~~~~~~~
When the final estimate (lambda)=(lambda) of the eigenvalue has
been found, the routine integrates the differential equation once
more with that value of (lambda), and with initial conditions
chosen so that the integral
b
/ 2 ddq
S= |y(x) ----------(x;(lambda))dx
/ dd(lambda)
a
is approximately one. When q(x;(lambda)) is of the form
(lambda)w(x)+q(x), which is the most common case, S represents
the square of the norm of y induced by the inner product
b
/
<f,g>= |f(x)g(x)w(x)dx,
/
a
with respect to which the eigenfunctions are mutually orthogonal.
This normalisation of y is only approximate, but experience shows
that S generally differs from unity by only one or two per cent.
During this final integration the REPORT routine supplied by the
user is called at each integration mesh point x. Sufficient
information is returned to permit the user to compute y(x) and
y'(x) for printing or plotting. D02KEF passes across to REPORT,
not y and y', but the Pruefer variables (beta), (phi) and (rho)
on which the numerical method is based. Their relationship to y
and y' is given by the equations
______ ( (rho)) ( (phi))
p(x)y'=\/(beta)exp( -----)cos( -----);
( 2 ) ( 2 )
1 ( (rho)) ( (phi))
y= --------exp( -----)sin( -----).
______ ( 2 ) ( 2 )
\/(beta)
For the theoretical basis of the numerical method to be valid,
the following conditions should hold on the coefficient
functions:
(a) p(x) must be non-zero and of one sign throughout the
interval (a,b); and,
ddq
(b) ---------- must be of one sign throughout (a,b) for all
dd(lambda)
relevant values of (lambda), and must not be identically
zero as x varies, for any (lambda).
Points of discontinuity in the functions p and q or their
derivatives are allowed, and should be included as 'break-points'
in the array XPOINT.
Parameters
d02kef
Required Input Arguments:
xpoint (:) real
coeffn function (User-Supplied)
bdyval function (User-Supplied)
k integer
tol real
elam real
monit function (User-Supplied)
report function (User-Supplied)
Optional Input Arguments: <Default>
hmax (2,:) real zeros(2,length(xpoint))
delam real 0.25*max(1,abs(elam))
match integer 0
maxfun integer 0
maxit integer 0
ifail integer -1
Output Arguments:
elam real
hmax (2,:) real
delam real
match integer
maxit integer
ifail integer